/*! \file */
/*
* Copyright (c) 2008-2011, Thomas Handorf
* 
* This file is part of libScopes.
* libScopes is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
* 
* libScopes is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU Lesser General Public License
* along with libScopes.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "sbml.h"

LIBSBML_CPP_NAMESPACE_USE

/*!
\brief import network from sbml file
\param net the network object
\param filename the name of the sbml file
adds all reactions to the provided net object. returns a map which links sbml species / reactions ids to compound / reaction indexes
*/

nodemap importSBML(Net& net, const char* filename, bool autoCat){
   SBMLDocument* document;
   SBMLReader reader;
  
   document = reader.readSBML(filename);
   document->printErrors(std::cerr);
   Model* model = document->getModel();
   
   map<string,unsigned int> nodes;
   unsigned int i;
   for(i=0; i < model->getNumSpecies(); i++){
      Species* s=model->getSpecies(i);
      nodes[s->getId()]=i;
      net.setName(i,s->getId());
   }
   VS subs;
   VS prds;
   VF subsst;
   VF prdsst;
   VS inh;
   VS cat;
   string id;
   for(i=0; i < model->getNumReactions(); i++){
      subs.clear();
      prds.clear();
      subsst.clear();
      prdsst.clear();
      inh.clear();
      cat.clear();
      Reaction* r=model->getReaction(i);
      id=r->getId();
      if (id.length()){
         if (nodes.count(id)){
            cout << "node id " << id <<" already exists"<< endl;
            abort();
         }
         nodes[id]=i;
      }
      unsigned int j,rn=r->getNumReactants();
      for(j=0; j < rn; j++) {
         SpeciesReference* rt = r->getReactant(j);
         double st=fabs(rt->getStoichiometry());
         subs.push_back(nodes[rt->getSpecies()]);
         subsst.push_back(st);
      }
      rn=r->getNumProducts();
      for(j=0; j < rn; j++) {
         SpeciesReference* rt = r->getProduct(j);
         double st=fabs(rt->getStoichiometry());
         prds.push_back(nodes[rt->getSpecies()]);
         prdsst.push_back(st);
      }
      rn=r->getNumModifiers();
      for(j=0; j < rn; j++) {
         ModifierSpeciesReference* rt = r->getModifier(j);
         int sbo=0;
         if (rt->isSetSBOTerm()) sbo=rt->getSBOTerm();
         if (sbo==13 || sbo==461 || sbo==460 || sbo==535 || sbo==533 || sbo==534 || ((sbo==0 || sbo==19) && autoCat)) cat.push_back(nodes[rt->getSpecies()]);  // catalysts and necessary stimulation
         if (sbo==20 || sbo==206 || sbo==207 || sbo==537 || sbo==536 ) inh.push_back(nodes[rt->getSpecies()]); // inhibitors
      }
      char rev=0;
      if (r->getReversible()) rev=1;
      net.replaceRea(i,&subs,&prds,rev,&subsst,&prdsst,&cat,&inh); // no inhibited substrate; no propagators
   }
   return nodes;
}
/*!
\brief import inhibitions definitions which are missing in the sbml files
\param net the network object
\param filename the name of the inhibition definitions file
\param nodes the nodemap generated by the sbml importer

retrives the inhibitions from the file. The inhibition definitions must use the same identifiers as the sbml file.
*/
void importInhibition(Net& net, const char* filename, nodemap nodes){
   cout << "importing inhibitions" << endl;
   FILE* f=fopen(filename,"r");
   if (!f){
      printf("could not find inhibition file %s\n",filename);
      abort();
   }
   char sr[255],si[255],sp[255];
   int ret=1;
   while (ret>0) ret=fscanf(f," #%[^\n]",sr); // remove comment lines
   ret=1;
   while (ret>0){ 
      ret=fscanf(f," %s -| %s + %s",sr,si,sp);
      if (ret>=2){ // inhibited substrate
         if (nodes.count(string(sr)) && nodes.count(string(si))){
            VS inhd;
            inhd.push_back(nodes[string(si)]);
            net.setInhibited(nodes[string(sr)],&inhd);
            cout << "found inhibited substrate '"<<si<<"' for reaction '"<<sr<<"'"<<endl;
         } else {
            cout << "at least one id has not been found in network'" << sr << "','"<< si << "'" <<endl;
            abort();
         }
      }
      if (ret==3){ //propagator (implies also inhibited substrate
         if (nodes.count(string(sp))){
            VS prp;
            prp.push_back(nodes[string(sp)]);
            net.setPropagators(nodes[string(sr)],&prp);
            cout <<"found inhibition propagator '"<<sp<<"' for reaction '"<<sr<<"'"<<endl;
         } else {
            cout << "id has not been found in network" << sp << endl;
            abort();
         }
         
      }
   }
}

/*!
\brief import family member definitions which are missing in the sbml files
\param net the network object
\param filename the name of the members definitions file
\param nodes the nodemap generated by the sbml importer

retrives the members from the file. The same identifiers as in the sbml file must be used.
*/
void importMembers(Net& net, const char* filename, nodemap nodes){
   cout << "importing memberships" << endl;
   FILE* f=fopen(filename,"r");
   if (!f){
      printf("could not find membership file %s\n",filename);
      abort();
   }
   char sf[255],sm[255];
   int ret=1;
   while (ret>0) ret=fscanf(f," #%[^\n]",sf); // remove comment lines
      ret=1;
   while (ret>0){ 
      ret=fscanf(f," %[^:]:",sf);
      if (ret>0){ // family 
         if (nodes.count(string(sf))){
            VS mem;
            while (ret>0){
               ret=fscanf(f," + %s",sm);
               if (ret>0){
                  if (nodes.count(string(sf))){
                     mem.push_back(nodes[string(sm)]);
                  } else {
                     cout << "id has not been found in network" << sm << endl;
                     abort();
                  }
               }
            }
            cout << "found members for '"<<sf<<"':";
            mem.print();
            net.setMembers(nodes[string(sf)],&mem);
         }
         else {
            cout << "id has not been found in network" << sf << endl;
            abort();
         }
         ret=1;
      }
   }
}


/*!
\brief import channel definitions which are missing in the sbml files
\param net the network object
\param filename the name of the channels definitions file
\param nodes the nodemap generated by the sbml importer

retrives the channels from the file. The same identifiers as the sbml file must be used for the compounds. Channels are just integers.
*/
void importChannels(Net& net, const char* filename, nodemap nodes){
   cout << "importing channels" << endl;
   FILE* f=fopen(filename,"r");
   if (!f){
      printf("could not find channel file %s\n",filename);
      abort();
   }
   char sf[255];
   int ret=1;
   while (ret>0) ret=fscanf(f," #%[^\n]",sf); // remove comment lines
      ret=1;
   while (ret>0){ 
      ret=fscanf(f," %[^:]:",sf);
      if (ret>0){ // family 
         if (nodes.count(string(sf))){
            VS mem;
            while (ret>0){
               int m;
               ret=fscanf(f," %d",&m);
               if (ret>0) mem.push_back((size_t) m);
            }
            cout << "found channels for '"<<sf<<"':";
            mem.print();
            net.setChannels(nodes[string(sf)],&mem);
         } else {
            cout << "id has not been found in network" << sf << endl;
            abort();
         }
         ret=1;
      }
   }
}
